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Abstract 

We study the nonlinear dynamics of the splitting of a doubly quantized vortex in a trapped 
condensate. The dynamics is studied in detail by solving the Gross-Pitaevskii equation. The main 
dynamical features are explained in terms of a nonlinear three-level system. We find an analytical 
solution for the characteristics of the dynamics. It is concluded that the time scale for the splitting 
is mainly determined by the instability of the linearized system, and nonlinear effects contribute 
logarithmically. 
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I. INTRODUCTION 



Quantization of fluid circulation is one of the most pictorial macroscopic manifestations 
of quantum mechanics. Lattices of singly quantized vortices have been imaged in supercon- 
ductors in magnetic fields ]| , liquid helium Q , and more recently in trapped Bose-Einstein 
condensates |3|, |4j. Vortices with higher quantum numbers than unity are energetically un- 
stable in many common situations, including an infinite, homogeneous s-wave superfluid, 
and the experimentally relevant case of a condensate contained in a parabolic potential 
, |6|, lZ|. In addition, in the latter case, multiply quantized vortices are found to be dynam- 



ically unstable in large areas of parameter space |8|, |9|, |l0|] . It was predicted that a doubly 
quantized vortex is unstable towards splitting into two vortices with unit quantum number, 
in accordance with the quantization of fluid circulation. 

These predictions were put to an experimental test in 2004, when a doubly quantized 
vortex, i.e., a vortex with quantum number 2, was imposed on a stationary condensate and 
the subsequent splitting was monitored Jj|. This experiment has been analyzed quantita- 



ivelyin Refs. 12j, |l3| using the time-dependent Gross-Pitaevskii equation [5], and in Refs. 

14] . \\\ by means of Bogoliubov theory. However, there remains to marry together these 
two approaches. In particular, Bogoliubov analysis gives information only about the linear 
(i. e., short-time) behavior of the unstable system, while solving the full Gross-Pitaevskii 
equation gives more detail than is necessary in order to understand the important features 
of the dynamics. 

Dynamics of vortices is a subject with a long history. It is well known that in a incom- 
pressible fluid the vortices move with the background fluid velocity [l^]. This is not so in 
a compressible fluid where the background density changes 161 ] . In general, vortex motion 



in a compressible fluid is complicated and cannot be separated from the dynamics of the 
system. The splitting of a doubly quantized vortex offers an opportunity to study the vortex 
dynamics in an extreme regime where the background velocity changes rapidly on the scale 
of the size of a vortex core. The splitting dynamics therefore offers insight into compressible 
fluid dynamics. In the study of the linear stability of doubly quantized vortices 15[], it was 
shown that the stability depends critically on the energy of the surface modes, and thus 
on global properties not associated with the vortex. The focus of this paper will be on the 
dynamics after the initial exponential growth of the vortex distance. Even though the ex- 
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periment of Ref. [llj was performed in an elongated three-dimensional geometry, this study 
is concerned with a two-dimensional system, in order to clearly bring out the structure of 
the problem. 

In this paper, we perform a systematic investigation of the long time behavior of the 
splitting of two vortices. The paper is organized as follows. In Sec. |TT]we discuss the equations 
governing the system. In Sec. II III we describe the numerical solution of the equations of 
motion. Section [IV] is devoted to a calculation of the nonlinear dynamics. The main features 
of the dynamics are captured in terms of a model that is solved analytically in Sec. |Vj Finally, 
in Sec. |VT]we summarize and conclude. Specifics of the analytical solution are given in the 
three appendices. 

II. SPLITTING OF A DOUBLY QUANTIZED VORTEX 

The system we study is a Bose-Einstein condensate of particles of mass m that is trapped 
in a cylindrically symmetric potential. At zero temperature in the dilute limit the gas 
is described by a condensate wavefunction ^(r, t) that obeys the Gross-Pitaevskii (GP) 
equation 

ih— = H V + U \V\ 2 V, (1) 

where 



Ht = -*Lr + vv, (2) 

and the trapping potential is assumed to be of the form 

V(r) = ^f(r 2 + X 2 z 2 ). (3) 

The inter-particle interactions are parametrized by an s-wave scattering length a, so that 
Uq = ATth 2 a/m. We immediately pass to trap units, where the unit of length is the oscillator 
length a osc = (h/mu) 1 ^ 2 and the unit of time is uj' 1 |5[. We assume the system to be two- 
dimensional (2D), which corresponds to the limit of a very tight trapping potential in the 
axial direction. The wavefunction in that direction is thus assumed to be in the ground state; 
on integrating out the z dependence one obtains the effective 2D interaction parameter 

C=— \Mz)\*dz = (4) 

«osc J a osc \J An 
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where 0o is the ground-state single-particle wave function in a one-dimensional harmonic 
potential. The resulting equation of motion for the condensate is 



-V 



-r 2 + C\V\ 2 
2 1 1 



(5) 



As a starting point for the study of the dynamics it is useful to repeat the linear 
stability analysis [l5]. The GP equation is expanded about a stationary solution 
^(r, t) = ^o(r) exp(— i/j,t) (which in the present case will be the doubly quantized vor- 
tex solution), where \x is the chemical potential of the system. The ansatz for the expansion 
is taken to be 



*(r,t) 



-i/it 



(6) 



where u n and v n are the quasiparticle amplitudes and u n the quasiparticle energies calculated 
from the Bogoliubov equations [5]. The small-amplitude excitations of the condensate are 
described by the eigenvectors and eigenvalues of the Bogoliubov equations, 



B 



where the linear operator B is defined by 



/ «n(r) ) 




' u n (r) \ 
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(7) 
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cm 2 



(8) 



If B has a complex eigenvalue, the system is dynamically unstable and the corresponding 
mode will grow exponentially. It is known that there exist intervals of the coupling constant 
C where the Bogoliubov equations possess a pair of complex eigenvalues. This behavior was 
thoroughly studied by the present authors in a previous paper 15] (cf. Q]). Figure CD shows 
the eigenvalue behavior as a function of C for the 2D case. An instability occurs when the 
energies of two Bogoliubov modes collide. In the present case the mode confined to the 
interior of the vortex, referred to as the core mode, mixes with surface modes of quadrupole 
symmetry. The core mode is seen in Fig. [T](a) as the line with positive slope that repeatedly 
merges with other lines representing the energies of quadrupole surface modes; each such 
collision creates an instability, so that successively higher instability regions correspond to 
increasing radial quantum number of the quadrupole mode. 
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FIG. 1: Energy levels in two dimensions for a condensate with an m = 2 vortex. The left panel 
shows the real parts and right panel the imaginary parts. All imaginary parts except at most one 
are zero at any point in this phase space. 

The Bogoliubov equations describe only the linear, i. e., small-amplitude, evolution of 
the condensate. In order to capture the full, nonlinear time development, in general one has 
to perform a numerical time integration of the time-dependent GP equation (CQ). However, 
the purpose of the present paper is to study to what extent a simplified approach, based on 
the solutions to the Bogoliubov equations, will suffice, and therefore we will in subsequent 
sections compare the full numerical results to the simplified model. To study the splitting 
dynamics of a doubly quantized vortex one needs to choose a perturbed doubly quantized 
vortex as initial condition. The doubly quantized vortex state is a stationary, rotationally 
symmetric solution of the GP equation (jHJ) of the form 

* 2 M) = / 2 (r)e* 2e , (9) 

where the real amplitude /2(f) obeys the equation 

1 /2(r) = /i/ 2 (r). (10) 



1/19 9 2 2 



2Wr^ + ^ + ^ ) + C|/2W|! 
As the initial condition for dynamical simulations one needs to add a perturbation to the 
doubly quantized vortex state. For defmiteness, we have chosen to use the ground-state 
harmonic oscillator wave function as a perturbation, but as long as the perturbation is 
small, its exact form does not matter for the long-time evolution, after it is exponentially 
inflated. 
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III. NUMERICAL METHOD 



The GP equation (jSJ) is solved using a Hermite mesh in both spatial directions, and the 
time evolution is done using a Strang splitting that makes use of the tensor product structure 



of the linear problem [18j. For a sufficiently large grid, in our case 100 x 100 points, we get 
conservation of angular momentum to one part in 10 6 . This symplectic method is nearly 
optimal for the problem at hand, which was crucial in order to be able to scan the parameter 
range and to analyze the subtle nonlinear dynamics in detail. 

The Bogoliubov equation is solved separately. Due to the cylindrical symmetry, it is 



reduced to a ID eigenvalue problem, whose solution was described in Ref. [15]. 

One of the most important quantities to be discussed in the following is the distance 
between two vortices in the numerical time evolution. To measure this distance, we first 
identify the spatial points r which fulfill the criteria |r| < 4 and |^(r)| < 0.15Max( 1 J r ); 
these are the points of low density in the interior of the system. Using these points we do a 
least-square fit to the form 

y(z)=A(z + z )(z-z )e^ 2 / 2 , (11) 

where z is short for z = x + iy. The fit is done with respect to the two constants A and 
zq. This fitting function describes two vortices placed symmetrically about the origin and 
is found to be an accurate approximation for the wavefunction at all times, in accordance 
with the expectation that the instability of a doubly quantized vortex results in the vortex 
splitting into two. The fit for the parameter zq gives the positions of the two vortices as zq 
and — zq, and the vortex distance is d = 2\zq\. A good fit is very difficult to achieve for small 
separations, since the least-square method minimization problem is then very shallow and 
small numerical errors in the wavefunction give significant contributions. A more reliable 
method to find the qualitative time evolution is to notice that in the weakly interacting 
limit, the squared length |z | 2 is approximately proportional to the population of the lowest 
harmonic-oscillator eigenstate (see [3], Eq. (23)). Therefore we project the wave function 
onto the eigenstates of the harmonic-oscillator potential, 

a n , m (t)= J ^(x,y,t)4> ntm (x,y)dxdy, (12) 

where <pn,m, is the eigenstate of the harmonic-oscillator potential with energy io n , m = 2n + 



6 



m| + 1, 




/ 



7r(n + m)\ 



L nm {r 2 ) (re 16 ) 




(13) 



and 




(14) 



is a generalized Laguerre polynomial. The population of an excited state is defined as 



The integral in Eq. f[T2"j) is calculated using the Gauss-Hermite quadrature rule associated 
with the Hermite mesh, which is exact in the limit of low energies. As we shall se, we find 
the amplitudes P n ,m, useful for understanding the dynamics of the problem. 

IV. TIME DEVELOPMENT OF VORTEX SPLITTING 



As known from previous studies [9|, |l5j, the dynamics of a perturbed doubly quantized 
vortex falls into one of two categories depending on the value of the coupling strength 
C . In some intervals the doubly quantized vortex is stable and in others it is unstable, as 
investigated in detail in Ref. [l5]. The real and imaginary parts of the Bogoliubov eigenvalues 
are presented in Fig. [TJ The regions where the vortex is linearly stable are not interesting 
from a dynamical perspective when small perturbations are considered. The condensate will 
just perform small periodic oscillations following the initial perturbation. Thus the domains 
of interest are the unstable regions. It turns out that these can roughly be divided into two: 
the first unstable region, and all the subsequent ones. 

We first consider the first unstable region, C G [0,37]. An example of the dynamics is 
given in Fig. [2j The depletion of the condensate, i. e. the m = 2 state, is very strong. It is 
seen that the sum of the populations in the m = 0, 2, and 4 states is less than 1 after some 
time, which means that there is a non-negligible population in states with m > 4. (Although 
the negligible population of states with odd m is here a consequence of the chosen initial 
conditions, we have checked that for more general initial conditions it is enforced by the 
dynamics, since only modes with even m become dynamically unstable.) The population 
in states with m > 4 is seeded by the large population in the m = 4 state, as will be 
clear below. Another feature which is worth noticing is that the vortex distance is highly 



Po,m (t) = |a , m (t)| • 



(15) 
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FIG. 2: (a) Time development of the vortex distance d, and (b) time development of the popula- 
tions P m o(t) of the harmonic-oscillator eigenstates, for a two-dimensional condensate with coupling 
parameter C = 20. In (b), the curves represent from the bottom up, m = 6, m = 4, m = 0, m = 2, 
and the sum of all four. The initial state was perturbed by means of a seeding of the harmonic 
oscillator ground state with an amplitude Poo = 0.001. 

correlated with the m = population, as anticipated in Sec. [IVl We take advantage of this 
near proportionality to find the time dependence of the vortex distance when the fitting 
method to find the vortex position described in Sec. [IV] fails. 

The time evolution proceeds in two stages. From the start the population of the m = 
state (which is the perturbation inserted by hand) and the m = 4 state grow exponentially 
while the condensate, the m = 2 state, is accordingly depleted. After the population of the 
m = and m = 4 states has become non- negligible, the populations of the two amplified 
states becomes asymmetric, due to population of higher-angular momentum states. The 
vortex distance and the population will start oscillating around finite values. Later we will 
see that the asymmetry and the population of higher- angular momentum eigenstates are 
crucial for the vortex distance to not oscillate back to zero. It is important to note that 
the asymmetry between the m = and m = 4 populations is not caused by the initial 
population chosen here, but is enforced by the dynamics. 

The dynamics in the higher unstable regions is different from that in the first. Figure [3] 
plots the population in different m states for C = 380, which is located in the third instability 
region (see Fig. [T]). Like in the first instability region, the time evolution of the unstable 
modes starts with an exponential growth. It achieves a maximum and start to oscillate. In 
contrast to the small-C case the oscillation is dominated by one frequency. Furthermore, 
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FIG. 3: (a) Square of the vortex distance d, and (b) total population in the m = (upper curve) 
and m = 4 (lower curve) states, for a condensate when the coupling parameter C = 380. The 
initial seeding of of the harmonic-oscillator ground state is Poo = 1 X 10 -5 . 

it is seen in Fig. [3] that the excited-state populations are very small at all times, and so is 
the depletion of the condensate. This is a general feature of the time evolution of higher 
instability regions, and it will enable us to make a simple model that captures the main 
features of the vortex dynamics and at the same time is analytically solvable (see Sec. IVj) . 
Finally, it is seen that the inter-vortex distance shows the same time dependence as the mode 
population Pq,o- The maximum distance between the two vortices is about d = 1 (in units 
of the oscillator length a osc as always), which is similar to that in the first unstable region, 
but contrary to that case, the diameter of the condensate is now much larger, meaning that 
the two vortices will stay well inside the condensate. The vortices rotate around each other 
and the distance between them oscillates in a non-sinusoidal way. 

From the discussion above, we may identify the most important characteristics of the 
dynamics as follows: (i) the exponential growth factor, (ii) the time until the first maximum 
is achieved, (iii) the maximum of the amplitude of the excited state, and (iv) the maximum 
inter-vortex distance. All of these features are functions of the nonlinear parameter C only. 
It is seen that items (i) and (ii) are closely related. The growth factor is given by the largest 
complex part of the Bogoliubov eigenvalues, while the time until first maximum must be 
inferred from numerical calculations; a comparison of these two is shown in Fig. HI The 
dashed line in Fig. H] is the result of the three-state model that will be described in Sec. 
M below. We see that in all instability regions the imaginary part of the mode frequency 
agrees well with the inverse of T max . 
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FIG. 4: Time for splitting of a doubly quantized vortex. Full lines represent the imaginary part of 
the complex Bogoliubov eigenvalue, asterisks represent the time T max taken for the vortex distance 
to achieve its first maximum according to the full GPE solution, and the dashed line is the same 
time scale in the three-state model, Eq. (|C15p . 

On the other hand, items (iii) and (iv), the maximum amplitude of the vortex distance 
and the maximum population of the unstable mode, show a quite surprising behavior. For 
the first instability region we see in Fig. [5] that both the maximum of the vortex distance 
and the maximum of the population are approximately independent of C in the unstable 
region. This is despite that the time to achieve this maximum varies strongly with C. 

For the higher unstable regions the behavior is very different. We present the result for 
the region C G [135, 205] in Fig. EJ The behavior of the maximum amplitude is particularly 
interesting since it has an almost linear increase from the start of the unstable region and 
reaches a maximum at the strong-coupling end of the unstable region, where it jumps dis- 
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FIG. 5: (a) The maximum radius of the motion of the vortices, and (b) the maximum occupation 
of the m = state, as a function of coupling parameter C with initial seeding -Po(O) = 0.01. The 
chosen range of C values lies in the first instability region. 
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FIG. 6: (a) The maximum radius of the motion of the vortices, and (b) the maximum probability 
of the m = state, Po, as a function of coupling parameter C with initial seeding -Po(O) = 0.001. 
The chosen range of C values lies in the second instability region. 

continuously to zero. This behavior will be explained in Sec. |V] where it is also shown that 
just outside this discontinuity a finite-amplitude perturbation can bring the system into the 
unstable region where the amplitude grows approximately to the maximum value achieved 
at the discontinuity. Finally we note in Fig. [6] that the maximal vortex distance d displays 
a similar behavior if one takes into account that Pq ~ d 2 . 
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V. THREE-MODE DYNAMICS OF VORTEX SPLITTING 



This section is devoted to extracting the main features of the dynamics of the splitting 
dynamics of the doubly quantized vortex that was studied numerically in the previous sec- 
tion. We put up a nonlinear model that can be solved analytically and which captures the 
main dynamics of the full system. The parameters of the model can be extracted from the 
the GP equation, for the most part analytically, and are all functions of C. This model will 
be particularly accurate for the higher unstable regions. 



A. First instability region 

It is instructive to first consider the dynamics in the first unstable region [19]. To find 
approximative solutions to the the GP equation, it is useful to start from the Lagrangian 
rom which the full GP equation can be derived if no further approximations are invoked 

r 1 d^* c 

L = i J dv-(**— - ¥— ) - (**# * + ^ |*| 4 ). (16) 

In the limit of small C, we know that the the dynamics mainly involves three states, namely 
the lowest-energy harmonic-oscillator eigenstates 4> n>m = 4>o,o, 4>o,2, and O ,4 ( see Eq-HS]). The 
m = 2 state represents the condensate, and m = and m = 4 the core and surface states 
respectively, which will be populated due to the instability. To investigate the dynamics of 
the vortex splitting in the space spanned by the three states, we expand the wave function 
as 

(r, t) = a o (t)0 o ,o(r) + a 2 (t)0 o , 2 + a 4 (t)0 o ,4, (17) 

so that a m is the amplitude of the state with m quanta of angular momentum in the in- 
direction. If we insert this into the Lagrangian we obtain 

/ * ■ ' * \ 

[ Ci' m Clrn &' 7n ( 3 j m ) 



= <E<« 

rn 

^ ^ ^m^rrflm ^ ^ ^ Cm,m ( 1 1m I ) 
m m 



I l 2 l I 2 

m<m' 

..2* , „*„2 * 



K(aoa 2 CL4 + CLQa 2 al) (18) 
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where C TO , m ' = C J \4>o,m\ 2 \4>o, m >\ 2 d 2 r, with m,m' = 0,2,4, and K = C J |0 o ,o| |0o,2| 2 |0o,4|d 2 r. 
The angular momentum conservation is automatically taken care of by the symmetries of the 
eigenfunctions. Using this expression we can write down the nonlinear equations of motion 



ia = (E + Co fi \a \ 2 + C 0: 2\a 2 \ 2 + Co^a^ao + Kajal (19) 
ia 2 = {E 2 + C 2fi \ao\ 2 + C 2 ^ 2 \a 2 \ 2 + C 2A \a 4 \ 2 )a 2 + 2KaQa* 2 a A (20) 
m 4 = (E + C^ao] 2 + C^ 2 \a 2 \ 2 + C iA \a 4 \ 2 )a4 + Ka^a 2 ,. (21) 

(22) 



Making a variable change a m — > a m = a m e iem<yt \ with a suitable choice of phases 9 m , the 
system of equations can be rewritten as 

ia = Ka 2 2 a\e im 
ia 2 = 2Ka Q a* 2 a i e~ i4 ' {t) 

ia 4 = Ka&le^, (23) 

where the phase is 

<pit) = j (E + C ,o|a | 2 + C , 2 \a 2 \ 2 + C 0A \a A \ 2 ) (24) 
+ (E A + C 4 , |ao| 2 + C 4i2 \a 2 \ 2 + C iA \a 4 \ 2 ) (25) 
- 2(E 2 + C 2 , |a | 2 + C %2 \a 2 \ 2 + C 2A \a A \ 2 )dt. (26) 

It can be seen from Eq. (1231 . that if initially ao(0) = 04(0), then it holds that ao(0) = 04(0) 
at all times. 

This equation is expected to give the correct dynamics for small C, where the dynamics 
is dominated by the lowest-energy single-particle eigenfunctions. For larger C values, the 
structure of the equation is still expected to be the same. As we will see below, the m = 2 
wavefunction should then be replaced by the condensate wavefunction and the m = and 
m = 4 states with the two Bogoliubov states associated with core excitations and surface 
excitations, respectively. 

In Fig. [7] we see that as long as the population is concentrated to the three states used 
in the truncation, the evolution of the truncated equation is identical to the full solution 
shown in Fig. [2j The main cause for discrepancies is that the m = 6 state starts to become 
populated, which causes a relative depletion of the m = 4 state. This results in an asymmetry 
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FIG. 7: The time dependence of the amplitude of the lowest m states in the truncated three-state 
model in equation IHol compared to the full GPE simulation. The interaction strength is C = 5 
and the initial seeding Pq = 0.1. From top to bottom, the full lines represent the amplitudes P2, 
Po, and P4 in the truncated model and the dotted lines represent the corresponding amplitudes for 
the full GPE simulation. 

between the populations of the m = 4 and m = states, and also implies that the amplitude 
ao (and with it the vortex distance) will not return back to zero, as it does in the truncated 
model. The terms in the Lagrangian (fT6l) causing this conversion are of the form G^a^a-oCte 
and a^a^ag, and thus they are proportional to three powers of the excited-state populations. 
In situations where the population of all higher states is small, the population of m = 6 and 
higher states is expected to be much slower, and this is also seen in Fig. [3j thus a three-state 
model should be more accurate for higher instability regions than for the first. 
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The truncated system has six degrees of freedom, corresponding to the real and imag- 
inary parts of the three amplitudes a m , but it has to conserve energy, norm and angular 
momentum, which leaves three degrees of freedom. In addition, the relative phase of the 
coefficients do and 04 will according to Eq. fliZ5|) stay constant; thus the system has only two 
degrees of freedom, which makes it integrable, so that the solution is periodic. 

The temporal dynamics depends on the initial state; since the initial increase of the 
excited-state population is exponential, it is expected that the time taken to achieve the first 
maximum is proportional to the logarithm of the initial population of the excited states. It 
is checked numerically that this logarithmic behavior is in fact very accurate even for the 
full nonlinear evolution. 

B. Higher instability regions 

For the dynamics in the higher unstable regions we have to modify the three-state model 
so that it takes into account the energy of the condensate and the coupling dependence of 
the quasiparticle energies. The structure of the model should also be such that it conserves 
angular momentum, quasiparticles and energy. The ansatz is therefore written 



tt(r, t) = e iftt [a 2 (t)* 2 (r) + a (t)u (r) + a*{t)v*{r) + a 4 (t)u 4 (r) + a^K(r)] , (27) 



where ^(r) is the condensate wavefunction with a doubly quantized vortex, and u and 
Vq are the Bogoliubov amplitudes for the core mode with m = 0. Finally, w 4 and V4 are 
the Bogoliubov amplitudes for a selected quadrupole mode, which is expected to become 



unstable when it mixes with the core mode. In Ref. [15] it was found that an instability occurs 
when the energy of the (n,m) = (0,0) Bogoliubov mode, the core mode, becomes nearly 
degenerate with a quadrupole mode with quantum numbers (n,m) = (n, m); the recurring 
instability regions arise from the crossings with quadrupole modes with successively higher n 
values; this can be seen in Fig. [Q All the functions in Eq. (j27j) are assumed to be calculated 
from Eqs. (15|71) at some fixed coupling strength C outside of any instability region; their 
energies are then to be extrapolated into the instability region. 

The calculations are carried out in Appendix [A] As already noted, the energies of the 
two Bogoliubov modes are nearly degenerate, and are assumed to coincide at a coupling Co. 
Furthermore, since the core mode is concentrated to the interior of the vortex, its energy 
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varies much more rapidly with coupling strength C than that of the quadrupole mode, so 
that only the C dependence of the former needs to be taken into account. Again, this is 
seen in Fig. [TJ Moreover, the same confinement also leads to a self-interaction of the core 
mode; corresponding terms for the other modes are small in comparison. Putting all this 
together results in the coupled equations 

ia = (2dE — 2I \a \ 2 )a + Ka\a\ 
ia 2 = 2i^a a 2 a 4 

ia 4 = Ka* a\. (28) 

Here, dE = (u — cg> 4 )/2 is half the C-dependent energy difference between the Bogoliubov 
eigenenergies; at the resonant coupling strength Cq we have 2dE{Co) = ujq(Cq) — ui^ = 0, so 
we may write 

dE =1^(0 -C„). (29) 

Inserting the 2D Thomas- Fermi approximation jlj] fi(C) = (C/7r) 1//2 , and using the expres- 
sion for the core mode energy [15], uo = 0.42/i, we obtain 

dE = 0A2—}==(C-C ). (30) 

The term J in Eq. (|2"8"|) represents the nonlinear self-interaction of the core mode, 

J = -C / dr\v \ 4 (r) w ° , , (31) 
where the last equality was carried out in App. [A] With Thomas-Fermi estimates for the 



core mode frequency 



^ (C) = 0.42J^, (32) 



and the quadrupole mode frequency for the n'th radially excited state [2( 

uo A = V2n 2 + 6n + 2, (33) 

the resonant coupling Cq was obtained in Ref. [jjj] as 

C = -^(2n 2 + Qn + 2), (34) 

where each value of n corresponds to an instability region. Finally, the constant K represents 
the integral that couples the three modes; it is found that any attempt to approximate this 
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term analytically is extremely sensitive to small variations in the variational parameters, so 
K has to be determined numerically. This can be done by noting (as will be shown in a 
moment) that the constant is fact equal to the maximum of the imaginary parts of the mode 
frequencies over the instability interval; numerically it is seen to be close to K m 0.15 for all 
instability regions. We note that Iq is at least an order of magnitude larger than K when C 
is of order 100 or more; this inequality will be taken advantage of in the calculations. Also 
note that whereas K and Iq are positive as long as the interactions are repulsive, the sign 
of dE depends on C — C . 

To see that K is related to the maximum imaginary part of the mode frequency, linearize 
Eq. (J2SJ) by removing the term proportional to J and put \a 2 \ = 1; the resulting oscillating 
solution for the amplitudes a and a 4 has a frequency 



in accordance with Bogoliubov theory. From this we conclude that the mode is unstable 
when \K/dE\ < 1, and that K is indeed the maximum imaginary part of the frequency. 

In Appendix [B] it is shown how the system of equations fl28|) leads to the differential 
equation for the core mode population p = |ao| 2 , 



p 2 = - [2(dE - K)p - (Jo - AK)p 2 - E ] [2{dE + K)p - (J + AK)p 2 - E ] = /(p), (36) 



This is an elliptic integral since f(p) is a polynomial of degree 4. The solution for p(t) is 
therefore given as an inverse of this elliptic integral. To understand the dynamics of the 
system we look at the zeros of f(p). The solution will oscillate between the two positive 
roots of f(p), since they correspond to p = 0. In the limit E <C (K, dE) <^ I [which holds 
according to the discussion below Eq. fl34l) ]. the roots can be written 



^lin — 



dE ± VdE 2 - K 2 



(35) 



where the constant Eq is the total energy. A formal solution is 




(37) 



Po = 



2(K + dE) 

dE + K 
2 ? 



2 ' 



Puiax 



-Po- 



(38) 



In a typical experimental situation, p is the initial value. 
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In App. Oit is found that the asymptotic expansion for the time to the first maximum, 
under the inequalities stated above, is given by 

1 



T 

k' 2 



=ln (16\ 

WK 2 - dE 2 n \k a J 
1 I 2 



(39) 



4 (l_(^)2)2(^ + rf £)2J 

The time scale is set by the imaginary part of the eigenvalue of the linearized problem, 
y/K 2 - dE 2 , as long as we have \dE/K\ < 1, as discussed in connection with Eq. (j3~5j) . The 
contribution of the initial population is only logarithmic. The nonlinearity described by 
the constant 1$ also contributes a logarithmic term. We conclude that the splitting time is 
mainly predicted by the linear Bogoliubov theory, and the nonlinear dynamics contributes 
only weakly. 

To further understand the dynamics is is useful to examine the Hamiltonian associated 
with these equations of motion. It is given by Eq. (1B12I) as 



H(ip d ,p) = 2K(- -p)pcosip d + (2dE - I p)p, 



(40) 



where tj) d is twice the phase difference between the two modes a 2 and a , = 2^2 — 2i>o- 
Again, we can according to Eqs. ( l3"0ll3~Ti) take the physically motivated limit K/Iq <^ 1. 
In Fig. [8] we see a contour plot of the Hamiltonian. Fig. E](a) shows the case dE < —K, 



a) 0.2 






FIG. 8: Contour plots of the Hamiltonian <[B12|) for dE = -1.5 (a), dE = (b), and dE = 1.5 (c). 
The other parameters are chosen as Iq = 10 and K = 1. The thick line is the separatrix, dividing 
phase space into running and oscillating solutions. 



where p = is a global maximum and the Hamiltonian is strictly convex, i.e., p = is 
unconditionally stable. In the linearly unstable regime, K > \dE\ [Fig.|S(b)], the point p = 
is a saddle point and the Hamiltonian has a global maximum for p = (dE + K)/Iq, ipd = 0. 
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The thick line in Fig. [8] is the separatrix, separating the solutions where ipd oscillates from 
the running ones. Finally, for dE > K, as shown in (c), p = is a local minimum and the 
Hamiltonian has a saddle point at p = (dE — K)/I and ipj = it. In this case the solution 
is stable when the initial conditions are sufficiently near p = 0; else it may start to oscillate 
around the maximum. 

The zeros of the polynomial f(p) defined in Eq. (|36|) correspond to points where the 
tangent of a contour line is horizontal. We observe that the contours are of two different 
types depending on whether they are closed lines, that do not wind about the origin, or 
whether they wind around phase space and connect at ipd = ±7r. If the initial condition is 
purely imaginary, ipj = tt, then the solution will always lie on a curve that winds around 
phase space. In the case ^ = the solution can lie on any level curve depending on the 
initial condition. Consider the case where dE > K, i.e., C is above the unstable region. 
Then a small initial value of p will yield a solution that lies in the stable region, i.e. the 
solution will circle around the local maximum. However, if the initial value of p is increased, 
the system enters a trajectory that winds around the minimum and p starts to oscillate. 
This is the reason for the finite-amplitude instabilities above the upper limit of the unstable 
region that were observed numerically in Sec. IIVI 

In Figs. I91TT01 and Fig. H] we compare the three-state model calculation with the full time 
integration. The overlap p is in the numerical calculation defined as the overlap integral 




t 



FIG. 9: Time development of the population p of the core mode. Full line represents the full 
numerical time integration, while the dashed line is obtained from the three-state approximation. 
The coupling strength is chosen as C = 380. 

between the condensate wave function and the core mode, analogously to Eq. f|T2|) . but with 



19 




_2 -1.5 -1 -0.5 0.5 1 1.5 2 

dE/K 

FIG. 10: Maximum amplitude po of the core mode as a function of energy difference divided by 
mode coupling, dE/K. The full line is the analytical solution of the three-state model, and the 
symbols represent data from the GPE solution in successive instability regions; circles for the 
second instability region, crosses for the third, and dots for the fourth. 

the numerically calculated core mode vq substituted for the single-particle eigenfunction o ,o- 
Fig. [10] collapses numerical data from the solution of the GP equation in several instability 
regions onto the same graph, by for each instability region mapping the coupling parameter 
C onto the parameters dE and K. As seen, the agreement with the three-state model for 
the time T in Fig. H]is excellent, but the magnitude of the amplitude p, seen in Fig. [101 is 
more sensitive to the exact parameter values, which may explain the discrepancy by a factor 
of order unity. Again, p does not oscillate back to zero in the full solution, but it does in the 
truncated model. Clearly, the population of higher states makes the dynamics nonperiodic, 
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so that the two vortices stay apart after they have separated. 



VI. CONCLUSIONS 

This study is concerned with the compressible vortex dynamics in a trapped Bose-Einstein 
condensate. The dynamics of the splitting of a doubly quantized vortex is studied in detail 
both in full numerical time integration, linear Bogoliubov analysis, and using a three-state 
model utilizing the Bogoliubov eigenstates. It is found that the simple three-state model 
captures many essential features of the dynamics. Moreover, it is seen that Bogoliubov anal- 
ysis is capable of determining the time scale for vortex splitting, while nonlinear processes 
only contribute logarithmically. 
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APPENDIX A: DERIVATION OF THREE-STATE MODEL 

We start from the ansatz 

(r, t) = e*** [a 2 (t)tf 2 (r) + a {t)u {r) + a* {t)v*{r) + a 4 {t)u A {r) + a^K(r)] , (Al) 

where u m , v m are the exact Bogoliubov amplitudes associated with the stationary condensate 
wave function ty 2 computed for a nonlinearity parameter C. C is assumed to lie outside of 
all instability regions, and the energies will be extrapolated into them. The dimensionless 
units were discussed in connection with Eq. @. We assume that the core state can be 
approximated well as a pure hole state. This assumption amounts to putting uq = 0. The 
functions thus fulfill the equations 

J d r %(r)(H + C\y 2 (r)\ 2 )y 2 (r) = fi, 
J drv*(r)(H + 2C\V 2 (r)\ 2 - n)v (r) = -u (C), 
J dvul(r)[(H + 2C|M> 2 (r)| 2 - ^(r) + C*!(r)i; 4 (r)] - 

J drvt(r)[(H + 2C|^ 2 (r)| 2 - ^(r) + C(^)V(r)] = u^C). (A2) 
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With the chosen sign conventions, uoq and uo 4 are both positive. In the absence of insta- 
bilities, a 4 is expected to oscillate as exp(— iu^t), and a Q oscillates as exp(— iuot). Taking 
nonlinearities into account, there will of course be corrections to the time dependence. Also 
note the dependencies on the azimuthal angle 9: ^2 exp(2i0), w 4 oc exp(4id), but t> and 
V4 do not depend on 9. 

On inserting the ansatz (1A1|) into the Lagrangian (fl6|) . it separates into five parts, 

L = L k + L + L r + L 2 + L 4 , (A3) 

where Lk contains the time derivatives, 

Lk = -(aQdo + 0^2 + — c.c.) 

+~(aod 4 J drv (r)t>* (r) - c.c). (A4) 

The term L contains the terms to which the eigenvalue equations flA2|) can be applied, 



L = -\a 2 \ 2 J dv[^2(r)\ 2 -^2(r)(H + C\^ 2 (r)\ 2 )^2(r)] 
-\a \ 2 J drv (r)(H + 2C|^ 2 (r)| 2 - fi)v*(r) 
-|a 4 | 2 f dr [ul(r)(H + 2C|^ 2 (r)| 2 - /i)u 4 (r) 

+v 4 (r)(H + 2C\y 2 (r)\ 2 - ^K(r)] , (A5) 



whereas L r collects the "rest terms" obtained from Lo because of the depletion of the con- 
densate, 



-C\a 2 \ 2 C-^f-l) J dr\^2{r)\' 
-2C|a | 2 (|a 2 | 2 -l) J rfr|M/ 2 (r)| 2 |i;o(r) 
-2C|a 4 | 2 (|a 2 | 2 -l) / dr|^ 2 (r)| 2 |n 4 (r) 



-2C|a 4 | 2 (|a 2 | 2 -I) J dr|^ 2 (r)| 2 K(r)| 2 . (A6) 
The term L2 contains terms of second order in the excited-state amplitudes, 

L 2 = -C[a* al(a 2 ) 2 J dr^ 2 (r)\*(rK(r) + cc] 
-Ca al\a 2 \ 2 J dr\^ 2 (r)\ 2 v (r)v i (r) 

-C[|a 4 | 2 |a 2 | 2 J dY^ 2 {v) 2 u\{v)v A {v) + c.c], (A7) 
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-TtKI 4 / '/r|r„(r) 



C 

— \aA z I (ir(k 4 (r)| 4 + IvJr) 



and analogously, L 4 contains the fourth-order terms 

C 
2 

_C 
~ 2 

— C|a | 2 |a 4 | 2 y c/r|t> | 2 (|t> 4 | 2 + | w 4 1 2 ) 

-C|a 4 | 4 y drK(r)| 4 Mr)| 2 . (A8) 

All the terms that are expected to oscillate as exp[i(u;o + uj i )t] are discarded. Next, consider 
all terms that are proportional to the fourth power of the excited-state amplitudes. Note 
that the function vq is concentrated to the vortex core, i.e., a very small spatial region, while 
the other functions, ^2, Ua, and t> 4 , are much less localized. As a result, the term in the 
first line in Eq. flA8j) is expected to be much larger than all the other terms in L 4 , and those 
are therefore discarded. Furthermore, all the terms in L r are also of fourth order in the 
excited-state amplitudes and can be discarded. The resulting Lagrangian is 

L = -(oqCIo + a 4 a 4 + a* 2 (i2 — c.c.) 

— \a \ 2 uj (C) — |a 4 | 2 a; 4 (C) 

— C ^aoa^al) 2 y c/ru 4 t>o | ^2 1 2 + c.c. 

— ^\ao\*J dv\v (v)\ 4 . (A9) 

Defining the constants 

|2 



K = C I dru 4 v Q \^ 
1 



2 



-C v, 



■ '0 • 



dE= l -[u Q {C)-u i {C)} 1 (A10) 
and making a phase change we can write the Lagrangian on the final form 

L = — (agdo + a 4 ^4 + a 2^2 — c.c.) 

- 2dE\a \ 2 - I \a \ 4 - (Ka Q a±{a* 2 ) 2 + c.c.) . (All) 

Following Ref. [l^], one may produce Thomas- Fermi estimates for Jo and dE, assuming 
that the core mode experiences an effective potential 



V {r)=2Cm 2 = 2^ 2 (A12) 



2 
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where £ = 1/ yflji is the healing length, and b is a variational parameter; the choice b = 2\/& 
minimizes the condensate energy. The ground state of this potential is 

l '° w = & xp( "^- (A13) 

and the nonlinear parameter of our model becomes 

In = *-= = = , (A14) 



where in the last line we used the Thomas- Fermi result \i — (C/tt 



iV2 



APPENDIX B: SOLUTION OF THE COUPLED NONLINEAR SYSTEM 

The Lagrangian for the system was in App. [A] found to be 

L(a m ,a* m ) = ^2 2 ( a m°>m ~ a m a* m ) - [(2dE - J K| 2 )K| 2 + K ((a* 2 ) 2 a a 4 + ala* a* 4 )] . 

m 

(Bl) 

First write the amplitudes in the form dj = exp(iipi). The Lagrangian can then be written 
in the form 

L(r m , ,f> m ) = ^r 2 m - {2Kr 2 2 r Q r A cos[2^ 2 - (Vo + M + (2dE - I Q r 2 )r 2 } . (B2) 

m 

Defining the auxiliary variables 

N = r 2 + r\ + rj, ^ n = 2^ 2 - 2^ 4 , 
L = 2r 2 2 + 4r 4 2 , ^ = -(^4-^0), 

D = r 2 Q + K\ + r 2 A , = -2^2 + V^o + ^4, (B3) 
the Lagrangian can be rewritten once more as 



L(N, L, D, Vi, fa) = i>nN + frL + faD - v F(D) cos^ d + G(D), (B4) 

where 

F(D) = AK 2 (2N -2D) 2 (D - L/A)(-N + L/A + D), 

G(D) = [2dE - I (D - L/4)](D - L/4). (B5) 
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We see that N and L are conserved; they are in fact the norm and angular momentum, 
respectively, so we suppress them as arguments. In addition the energy function 

E = y/F{D) cos + G(D) (B6) 

is conserved. The Lagrange equations for D and ipd read 

D = -y / F(D) sin^d, 

4d = w ^/F\D) C o^ d + QpG(D). (B7) 

Using energy conservation and the square of the first line we obtain the ODE 

D 2 = F(D) — (E — G(D)) 2 . (B8) 

This is an elliptic ODE since the right hand side is a polynomial of degree 4. The solution is 
done by factorizing the polynomial on the right-hand side into two second-order polynomials. 

Suppose N — 1, L — 2, which is the case in the present physical problem. In this case 
F(D) is a square, 

F(D) = f{Df = 4K 2 (1 - Df{D - i) 2 , (B9) 
and the equation (IB8I) for D simplifies to 

D 2 = (f(D) + E- G{D)){f{D) — (E — G{D))). (BIO) 

Now rewrite the equation in terms of the original variable p = r 2 and obtain the final 
equation of motion, 

p 2 = P\{p) P2 {p) , where 
Px(p) = (-4K + I )p 2 + 2(K-dE)p + E, 

P 2 (p) = (-4K-I )p 2 + 2(K + dE)p-E. (Bll) 

Note that the roots of the polynomial P\(p)P2{p) may be either real or complex; two roots 
will become complex when (dE — K) 2 < (I — 4K)E. Since E is proportional to the initial 
population po [see Eq. (1381) ]. it can be assumed small and hence the complex roots appear 
only in a very small portion of phase space; this permits us to concentrate on the case with 
real roots only. Furthermore, as is shown in Sec. |Vj we can on physical grounds assume 
Iq K; this will simplify some expressions in the following. 
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It is useful to write these equations as a Hamiltonian system with ipd,P canonically 
conjugate variables. Starting from Eqs. flB7j) for ipj, D and shifting variables as above to 
ipdiV we obtain 

H(if> d ,p) = 2K(±- p)p cosVd + (2dE - I oP )p. (B12) 
APPENDIX C: ELLIPTIC INTEGRALS 

We now solve the differential equation for the core- mode amplitude p, Eq. (1B11I) . Rewrit- 
ing this as 

P = Io\/(P- Pmin) (P + Pl) (Pmax - P) (j>2 + P) , (CI) 

where p m ; n and p max are the smallest and largest positive roots of the polynomial, respec- 
tively, and —pi and — p<i are the other two, and taking the initial value for p to be at the 
minimum point, then we may write 

t = I f m d _E . (C2) 

Now define p x = (p min + pi)/2,dl = (p min -pi)/2, p b = Qo max + p 2 )/2, db = (p max -p 2 )/2, 
and p' = p — dl, to obtain the integral 

t = i r d £ fC3 ) 

" J Pl \{P' 2 - P 2 )((p b + db- dl) - p>)((p b -db- dl) + p')] 1/2 ' 1 ' 



2l|, 



The asymmetry between the largest zeros can be removed by invoking the substitution 
p. 514) 

/ ax + b . . 

P = C4 
cx + a 

The parameters are to be determined so that the transformation leaves the symmetric zeros 
of the integrand invariant but makes the other two symmetric in terms of the new variable 
x; the new zeros of the denominator are denoted by ±p s . The integral for t is now 

t = ± r dx (cs) 

h J Pl v/(x 2 -p 2 )(p 2 s -x 2 Y 
and assuming that all roots are real, as discussed in App. (Bj the solution is 

1 &'4 = ri- i ), (C6) 



^ \ Pi \P 
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where nd is a Jacobian elliptic function ( 22j, p. 596), and 



A 2 = (arf - cdf 

((api) 2 - c 2 )(c - (a(p b + dr)))(c + a(p b - dr))' 

with dr = db — dl. This gives the complete result 

x(t) = pind(^-t), 
., ax(t) + b 

The half period T of the function p(t) is given by the complete elliptic integral 

T = ^-JC(k> 2 ), £;' 2 = (-) 2 - (C9) 
IqPs Ps 

Expanding the parameters in powers of pi, which in our physical situation corresponds 

to assuming that the initial population is very small, yields 



,dr 2 dr 2 1 pf 



^ « (1 - (-)') " — TTTT^ + W). (CIO) 

Vb - 



The transformation is given by 

\ c?r dr(dr 2 — p 2 )^ 1 ~^ ^ ^ 
6 = cp 2 . (Oil) 

Finally, the transformed root is 

P S = + , f h , a , rf + 0(pf). (C12) 

Pfe Pb(p 2 b -dr 2 ) 

The asymptotic expression for the complete elliptic integral when its argument is small is 

*W~ W^Y (C13) 



2 V fc ' 2 , 

Wrapping up all of the above, we obtain the full time evolution p(t) from Eq. (108|) where 
we substitute 

. n , dE . r, 

A 1 " { ~K ] • 

P. = £(1 - (f )*), (CM, 
and po is the initial population; the time taken to attain the first maximum is given by 

1 mf^Y (ci5) 



4^/1 -(f) 2 ^ 
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where 



1 I 2 

C ' 2 w (l-(f )2)2 4(ZTW P "' (C16) 
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